Arcsine Distribution

The arcsine distribution is a continuous distribution on a finite interval whose density spikes near the endpoints.

It shows up in two classic ways:

1. As a special Beta distribution: (\text{Beta}(\tfrac12,\tfrac12)) on ([0,1]).

2. As the limit law in the arcsine laws for symmetric random walks / Brownian motion (e.g., the fraction of time spent above 0).

---

Learning goals

By the end, you should be able to:

• write down and interpret the pdf/cdf of the arcsine distribution

• derive its mean and variance using a clean trig substitution

• sample it efficiently with NumPy only (inverse CDF)

• use scipy.stats.arcsine for pdf, cdf, rvs, and fit

• recognize where it appears in statistics (Jeffreys prior) and stochastic processes (arcsine laws)

Prerequisites

• Calculus: change of variables, basic integrals

• Probability: pdf/cdf, expectation and variance

• Familiarity with the Beta distribution is helpful but not required

import numpy as np
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio

from scipy import special, stats

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(42)

1) Title & Classification

Distribution name: arcsine

Type: continuous

Support:

• Canonical form: (x \in [0,1])

• General (endpoint) form: (x \in [a,b]) with (a<b)

Parameter space:

• Endpoint form: ((a,b) \in \mathbb{R}^2) with (a<b)

• SciPy form (location/scale): loc = a, scale = b-a with scale > 0

Key identity:

[ X \sim \text{arcsine on }[0,1] \quad\Longleftrightarrow\quad X \sim \text{Beta}\left(\tfrac12,\tfrac12\right). ]

2) Intuition & Motivation

What it models

The arcsine distribution is a natural model for fractions and time proportions that tend to be near 0 or near 1 more often than “in the middle”.

Its pdf is U-shaped on ([0,1]):

• very high density near 0 and 1

• lowest density near 1/2

This makes it a good model when boundary behavior is common (e.g., “almost always” vs “almost never”).

Where it comes from (a memorable construction)

Let (\Theta \sim \mathrm{Unif}(0,\pi)). Then

[ Y = \cos\Theta \in [-1,1] ]

has density (f_Y(y)=\tfrac{1}{\pi\sqrt{1-y^2}}), the classic arcsine law on ([-1,1]).

If we map ([-1,1]) to ([0,1]) via (X = \tfrac{Y+1}{2}), we obtain

[ f_X(x) = \frac{1}{\pi\sqrt{x(1-x)}}, \qquad x\in(0,1). ]

This trig construction is also the cleanest route to the mean/variance derivations later.

Typical real-world use cases

• Stochastic processes (arcsine laws): for a symmetric random walk / Brownian motion, the fraction of time the process stays above 0 has an arcsine limit law.

• Bayesian statistics: (\text{Beta}(\tfrac12,\tfrac12)) is the Jeffreys prior for a Bernoulli/Binomial probability parameter (p). It is “noninformative” in an information-geometric sense and is heavier near 0 and 1 than the uniform prior.

• Generative modeling of probabilities: when you want random probabilities that are often extreme (near 0 or 1), arcsine is a simple choice.

Relations to other distributions

• Special case of the Beta family: (\text{arcsine} = \text{Beta}(\tfrac12,\tfrac12))

• Location/scale transform to any ([a,b])

• Closely related to the “arcsine on ([-1,1])” density (\tfrac{1}{\pi\sqrt{1-y^2}})

3) Formal Definition

We’ll use the endpoint parameterization ((a,b)) with (a<b).

PDF

[ f(x\mid a,b) = \frac{1}{\pi\sqrt{(x-a)(b-x)}}, \qquad x\in(a,b) ]

and (f(x\mid a,b)=0) for (x\notin[a,b]).

Notes:

• The pdf diverges as (x\to a) or (x\to b), but the divergence is integrable.

• In SciPy, this corresponds to stats.arcsine(loc=a, scale=b-a).

CDF

[ F(x\mid a,b) = \begin{cases} 0, & x\le a,\ \frac{2}{\pi}\arcsin!\Big(\sqrt{\frac{x-a}{b-a}}\Big), & a<x<b,\ 1, & x\ge b. \end{cases} ]

In the canonical ([0,1]) case:

[ f(x)=\frac{1}{\pi\sqrt{x(1-x)}},\quad F(x)=\frac{2}{\pi}\arcsin(\sqrt{x}). ]

def _check_ab(a: float, b: float) -> None:
    if not (np.isfinite(a) and np.isfinite(b) and a < b):
        raise ValueError("Require finite endpoints with a < b.")


def arcsine_pdf(x, a: float = 0.0, b: float = 1.0):
    '''Arcsine pdf on [a,b]. Vectorized over x.'''
    _check_ab(a, b)
    x = np.asarray(x, dtype=float)

    out = np.zeros_like(x, dtype=float)
    interior = (x > a) & (x < b)
    boundary = (x == a) | (x == b)

    out[boundary] = np.inf
    out[interior] = 1.0 / (np.pi * np.sqrt((x[interior] - a) * (b - x[interior])))
    return out


def arcsine_logpdf(x, a: float = 0.0, b: float = 1.0):
    '''Numerically friendlier log-pdf (still +inf at boundaries).'''
    _check_ab(a, b)
    x = np.asarray(x, dtype=float)

    out = np.full_like(x, -np.inf, dtype=float)
    interior = (x > a) & (x < b)
    boundary = (x == a) | (x == b)

    out[boundary] = np.inf
    out[interior] = (
        -np.log(np.pi)
        - 0.5 * np.log(x[interior] - a)
        - 0.5 * np.log(b - x[interior])
    )
    return out


def arcsine_cdf(x, a: float = 0.0, b: float = 1.0):
    '''Arcsine CDF on [a,b].'''
    _check_ab(a, b)
    x = np.asarray(x, dtype=float)

    z = (x - a) / (b - a)
    z = np.clip(z, 0.0, 1.0)  # protects sqrt/arcsin from tiny floating-point drift

    out = (2.0 / np.pi) * np.arcsin(np.sqrt(z))
    out = np.where(x <= a, 0.0, out)
    out = np.where(x >= b, 1.0, out)
    return out


def arcsine_ppf(u, a: float = 0.0, b: float = 1.0):
    '''Inverse CDF (percent point function).'''
    _check_ab(a, b)
    u = np.asarray(u, dtype=float)
    if np.any((u < 0.0) | (u > 1.0)):
        raise ValueError("u must be in [0,1].")
    return a + (b - a) * np.sin(0.5 * np.pi * u) ** 2


def arcsine_rvs(size=None, a: float = 0.0, b: float = 1.0, rng: np.random.Generator | None = None):
    '''NumPy-only sampling via inverse transform.'''
    _check_ab(a, b)
    if rng is None:
        rng = np.random.default_rng()
    u = rng.random(size=size)
    return arcsine_ppf(u, a=a, b=b)


# Quick sanity check against SciPy on [0,1]
x_grid = np.linspace(1e-6, 1 - 1e-6, 5)
np.c_[x_grid, arcsine_pdf(x_grid), stats.arcsine.pdf(x_grid)]

array([[  0.000001, 318.310045, 318.310045],
       [  0.25    ,   0.735105,   0.735105],
       [  0.5     ,   0.63662 ,   0.63662 ],
       [  0.75    ,   0.735105,   0.735105],
       [  0.999999, 318.310045, 318.310045]])

4) Moments & Properties

Let (X \sim \mathrm{Arcsine}(a,b)).

Mean and variance

[ \mathbb{E}[X] = \frac{a+b}{2}, \qquad \mathrm{Var}(X) = \frac{(b-a)^2}{8}. ]

Because the distribution is symmetric around (\tfrac{a+b}{2}), the skewness is 0.

Kurtosis

For the canonical ([0,1]) case (equivalently (\text{Beta}(\tfrac12,\tfrac12))):

• skewness = 0

• excess kurtosis = (-\tfrac{3}{2}) (so kurtosis = (3 - \tfrac{3}{2} = \tfrac{3}{2}))

These values stay the same under location/scale transforms.

MGF and characteristic function

It’s convenient to express these using Bessel functions. Define:

• (I_0): modified Bessel function of the first kind (order 0)

• (J_0): Bessel function of the first kind (order 0)

Then

[ M_X(t) = \mathbb{E}[e^{tX}] = \exp\Big( t,\tfrac{a+b}{2}\Big), I_0\Big( t,\tfrac{b-a}{2}\Big), \qquad t\in\mathbb{R} ]

[ \varphi_X(t) = \mathbb{E}[e^{itX}] = \exp\Big( it,\tfrac{a+b}{2}\Big), J_0\Big( t,\tfrac{b-a}{2}\Big), \qquad t\in\mathbb{R}. ]

Entropy (differential)

For the canonical ([0,1]) distribution:

[ H(X) = \log\Big(\frac{\pi}{4}\Big). ]

On ([a,b]), scaling adds (\log(b-a)):

[ H(X) = \log(b-a) + \log\Big(\frac{\pi}{4}\Big) = \log\Big( (b-a),\frac{\pi}{4}\Big). ]

def arcsine_mean(a: float = 0.0, b: float = 1.0) -> float:
    _check_ab(a, b)
    return 0.5 * (a + b)


def arcsine_var(a: float = 0.0, b: float = 1.0) -> float:
    _check_ab(a, b)
    return (b - a) ** 2 / 8.0


def arcsine_mgf(t, a: float = 0.0, b: float = 1.0):
    '''MGF using modified Bessel I0.'''
    _check_ab(a, b)
    t = np.asarray(t, dtype=float)
    return np.exp(t * 0.5 * (a + b)) * special.i0(t * 0.5 * (b - a))


def arcsine_cf(t, a: float = 0.0, b: float = 1.0):
    '''Characteristic function using Bessel J0.'''
    _check_ab(a, b)
    t = np.asarray(t, dtype=float)
    return np.exp(1j * t * 0.5 * (a + b)) * special.j0(t * 0.5 * (b - a))


def arcsine_entropy(a: float = 0.0, b: float = 1.0) -> float:
    _check_ab(a, b)
    return float(np.log((b - a) * np.pi / 4.0))


# Monte Carlo check on [0,1]
n_mc = 200_000
x_mc = arcsine_rvs(n_mc, rng=rng)

print("Monte Carlo mean:", x_mc.mean(), "| theory:", arcsine_mean())
print("Monte Carlo var :", x_mc.var(), "| theory:", arcsine_var())
print("Entropy theory  :", arcsine_entropy())

Monte Carlo mean: 0.4999078439779247 | theory: 0.5
Monte Carlo var : 0.12498062340322307 | theory: 0.125
Entropy theory  : -0.2415644752704905

5) Parameter Interpretation

The parameters (a) and (b) are endpoints of the support.

• (a) shifts the distribution left/right

• (b-a) scales (stretches) the interval

Crucially, the shape in normalized coordinates does not change.

If (X \sim \mathrm{Arcsine}(a,b)) and

[ Z = \frac{X-a}{b-a}, ]

then (Z \sim \mathrm{Arcsine}(0,1)) (equivalently (\text{Beta}(\tfrac12,\tfrac12))).

So the family is “rigid”: parameters only translate and rescale.

intervals = [(0.0, 1.0), (-1.0, 1.0), (2.0, 5.0)]

fig = go.Figure()

for a, b in intervals:
    x = np.linspace(a + 1e-4 * (b - a), b - 1e-4 * (b - a), 800)
    fig.add_trace(
        go.Scatter(
            x=x,
            y=arcsine_pdf(x, a=a, b=b),
            name=f"pdf on [{a:g}, {b:g}]",
        )
    )

fig.update_layout(
    title="Arcsine pdf for different endpoints (note the endpoint spikes)",
    xaxis_title="x",
    yaxis_title="density",
)
fig.show()

6) Derivations

6.1 Expectation and variance (trig substitution)

A standard trick is to parameterize (x\in[a,b]) by an angle.

Let

[ x(\theta) = \frac{a+b}{2} + \frac{b-a}{2}\cos\theta,\qquad \theta\in(0,\pi). ]

Then

[ (x-a)(b-x) = \Big(\tfrac{b-a}{2}\Big)^2 \sin^2\theta, \qquad \mathrm{d}x = -\tfrac{b-a}{2}\sin\theta,\mathrm{d}\theta. ]

Plugging into (f(x\mid a,b),\mathrm{d}x):

[ \frac{1}{\pi\sqrt{(x-a)(b-x)}},\mathrm{d}x = \frac{1}{\pi,(\tfrac{b-a}{2})\sin\theta},\Big(-\tfrac{b-a}{2}\sin\theta,\mathrm{d}\theta\Big) = -\frac{1}{\pi},\mathrm{d}\theta. ]

Flipping integration limits turns the minus sign into a plus, so effectively

[ f(x),\mathrm{d}x = \frac{1}{\pi},\mathrm{d}\theta, \quad \theta\sim\mathrm{Unif}(0,\pi). ]

Mean.

[ \mathbb{E}[X] = \frac{1}{\pi}\int_0^\pi x(\theta),\mathrm{d}\theta = \frac{1}{\pi}\int_0^\pi \Big(\tfrac{a+b}{2} + \tfrac{b-a}{2}\cos\theta\Big),\mathrm{d}\theta = \frac{a+b}{2} ]

because (\int_0^\pi \cos\theta,\mathrm{d}\theta = 0).

Variance. Write (\mu=\tfrac{a+b}{2}). Then

[ X-\mu = \tfrac{b-a}{2}\cos\theta, \qquad (X-\mu)^2 = \Big(\tfrac{b-a}{2}\Big)^2 \cos^2\theta. ]

Therefore

[ \mathrm{Var}(X) = \mathbb{E}[(X-\mu)^2] = \Big(\tfrac{b-a}{2}\Big)^2 \cdot \frac{1}{\pi}\int_0^\pi \cos^2\theta,\mathrm{d}\theta = \Big(\tfrac{b-a}{2}\Big)^2 \cdot \frac{1}{2} = \frac{(b-a)^2}{8}. ]

6.2 Likelihood

For i.i.d. data (x_1,\dots,x_n) with all points in ((a,b)), the likelihood is

[ L(a,b) = \prod_{i=1}^n \frac{1}{\pi\sqrt{(x_i-a)(b-x_i)}}. ]

The log-likelihood (up to an additive constant) is

[ \ell(a,b) = -\frac{1}{2}\sum_{i=1}^n \log(x_i-a) - \frac{1}{2}\sum_{i=1}^n \log(b-x_i), \qquad \text{subject to } a<\min_i x_i,; b>\max_i x_i. ]

Important: because (\log(x_i-a)\to -\infty) as (a\uparrow \min x_i), the log-likelihood can grow without bound. So the unconstrained MLE for ((a,b)) does not exist (the likelihood is unbounded). Practical fitting methods therefore use constraints/regularization or alternative estimators.

7) Sampling & Simulation

Inverse transform sampling (NumPy-only)

Starting from the CDF for ([a,b]):

[ u = F(x\mid a,b) = \frac{2}{\pi}\arcsin!\Big(\sqrt{\frac{x-a}{b-a}}\Big) ]

Solve for (x):

[ \sqrt{\frac{x-a}{b-a}} = \sin\Big(\frac{\pi u}{2}\Big) \quad\Rightarrow\quad x = a + (b-a),\sin^2\Big(\frac{\pi u}{2}\Big). ]

Algorithm

1. Sample (u \sim \mathrm{Unif}(0,1))

2. Return (x = a + (b-a)\sin^2(\pi u/2))

This is exact, fast, and requires only NumPy.

# Sampling demo
a, b = 0.0, 1.0
x = arcsine_rvs(10, a=a, b=b, rng=rng)
x

array([0.35102 , 0.333278, 0.957626, 0.620713, 0.864426, 0.331615,
       0.148247, 0.814071, 0.448362, 0.69346 ])

8) Visualization

We’ll visualize:

• the pdf

• the cdf

• Monte Carlo samples vs the theoretical pdf

# PDF and CDF on [0,1]
eps = 1e-4
x = np.linspace(eps, 1 - eps, 1000)

fig_pdf = go.Figure()
fig_pdf.add_trace(go.Scatter(x=x, y=arcsine_pdf(x), name="arcsine pdf"))
fig_pdf.add_trace(go.Scatter(x=x, y=np.ones_like(x), name="uniform(0,1) pdf", line=dict(dash="dash")))
fig_pdf.update_layout(title="PDF on [0,1] (arcsine vs uniform)", xaxis_title="x", yaxis_title="density")
fig_pdf.show()

fig_cdf = go.Figure()
fig_cdf.add_trace(go.Scatter(x=x, y=arcsine_cdf(x), name="arcsine cdf"))
fig_cdf.add_trace(go.Scatter(x=x, y=x, name="uniform(0,1) cdf", line=dict(dash="dash")))
fig_cdf.update_layout(title="CDF on [0,1] (arcsine vs uniform)", xaxis_title="x", yaxis_title="F(x)")
fig_cdf.show()

# Monte Carlo samples vs pdf
n = 60_000
samples = arcsine_rvs(n, rng=rng)

fig = go.Figure()
fig.add_trace(
    go.Histogram(
        x=samples,
        nbinsx=80,
        histnorm="probability density",
        name="samples (hist)",
        opacity=0.6,
    )
)
fig.add_trace(go.Scatter(x=x, y=arcsine_pdf(x), name="theoretical pdf", line=dict(color="black")))
fig.update_layout(
    title="Monte Carlo histogram with theoretical pdf overlay",
    xaxis_title="x",
    yaxis_title="density",
    barmode="overlay",
)
fig.show()

print("sample mean/var:", samples.mean(), samples.var())

sample mean/var: 0.5012310506354244 0.12544470668500474

9) SciPy Integration

SciPy provides the distribution as scipy.stats.arcsine.

• stats.arcsine.pdf(x, loc=a, scale=b-a)

• stats.arcsine.cdf(x, loc=a, scale=b-a)

• stats.arcsine.rvs(size=..., loc=a, scale=b-a, random_state=...)

• stats.arcsine.fit(data) estimates loc and scale

We’ll also verify the identity with (\text{Beta}(\tfrac12,\tfrac12)).

# SciPy: pdf/cdf/rvs
x = np.linspace(1e-4, 1 - 1e-4, 1000)

pdf_scipy = stats.arcsine.pdf(x)
cdf_scipy = stats.arcsine.cdf(x)

# Identity: arcsine == Beta(1/2, 1/2)
pdf_beta = stats.beta(a=0.5, b=0.5).pdf(x)

print("max |pdf_scipy - pdf_beta|:", np.max(np.abs(pdf_scipy - pdf_beta)))

# Sampling
s_scipy = stats.arcsine.rvs(size=5, random_state=rng)
s_numpy = arcsine_rvs(5, rng=rng)
print("SciPy rvs:", s_scipy)
print("NumPy rvs:", s_numpy)

# Fitting loc/scale (note: likelihood is tricky near the endpoints)
data = arcsine_rvs(2_000, a=2.0, b=5.0, rng=rng)
loc_hat, scale_hat = stats.arcsine.fit(data)
print("fit loc, scale:", loc_hat, scale_hat)
print("true loc, scale:", 2.0, 3.0)

max |pdf_scipy - pdf_beta|: 3.552713678800501e-15
SciPy rvs: [0.425294 0.712055 0.864634 0.234072 0.38168 ]
NumPy rvs: [0.638791 0.463073 0.730777 0.226557 0.552068]
fit loc, scale: 2.0000041654452665 3.00059414578592
true loc, scale: 2.0 3.0

10) Statistical Use Cases

10.1 Hypothesis testing (arcsine law for random walks)

For a symmetric random walk (S_t = \sum_{i=1}^t \epsilon_i) with (\epsilon_i\in{-1,+1}) i.i.d., one version of the arcsine law says:

the fraction of time the walk is positive converges in distribution to an arcsine law.

We’ll simulate many random walks, compute the proportion of steps with (S_t>0), and compare to the arcsine distribution on ([0,1]).

n_paths = 4000
n_steps = 600

steps = rng.choice([-1, 1], size=(n_paths, n_steps))
paths = np.cumsum(steps, axis=1)

frac_positive = (paths > 0).mean(axis=1)

ks = stats.kstest(frac_positive, stats.arcsine.cdf)
print("KS statistic:", ks.statistic)
print("KS p-value   :", ks.pvalue)

x = np.linspace(1e-4, 1 - 1e-4, 800)

fig = go.Figure()
fig.add_trace(
    go.Histogram(
        x=frac_positive,
        nbinsx=60,
        histnorm="probability density",
        name="random-walk fractions",
        opacity=0.6,
    )
)
fig.add_trace(go.Scatter(x=x, y=stats.arcsine.pdf(x), name="arcsine pdf", line=dict(color="black")))
fig.update_layout(
    title="Fraction of time positive in a symmetric random walk (simulation)",
    xaxis_title="fraction of steps with S_t > 0",
    yaxis_title="density",
    barmode="overlay",
)
fig.show()

KS statistic: 0.02575
KS p-value   : 0.009764424160149908

10.2 Bayesian modeling (Jeffreys prior for Bernoulli/Binomial)

For a Bernoulli/Binomial success probability (p), the Jeffreys prior is

[ p \sim \mathrm{Beta}(\tfrac12, \tfrac12), ]

which is exactly the arcsine distribution on ([0,1]).

If we observe (k) successes in (n) trials, the posterior is

[ p \mid k \sim \mathrm{Beta}\Big(k+\tfrac12,; n-k+\tfrac12\Big). ]

We’ll compare Jeffreys’ prior to a uniform prior (\mathrm{Beta}(1,1)).

n, k = 20, 3

prior_jeffreys = stats.beta(0.5, 0.5)
prior_uniform = stats.beta(1.0, 1.0)

post_jeffreys = stats.beta(k + 0.5, n - k + 0.5)
post_uniform = stats.beta(k + 1.0, n - k + 1.0)

x = np.linspace(1e-4, 1 - 1e-4, 1000)

fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=prior_jeffreys.pdf(x), name="Jeffreys prior Beta(1/2,1/2)", line=dict(dash="dash")))
fig.add_trace(go.Scatter(x=x, y=prior_uniform.pdf(x), name="Uniform prior Beta(1,1)", line=dict(dash="dash")))
fig.add_trace(go.Scatter(x=x, y=post_jeffreys.pdf(x), name=f"Posterior (Jeffreys), k={k}, n={n}"))
fig.add_trace(go.Scatter(x=x, y=post_uniform.pdf(x), name=f"Posterior (uniform), k={k}, n={n}"))

ci_low, ci_high = post_jeffreys.ppf([0.025, 0.975])
fig.add_vline(x=ci_low, line=dict(color="black", dash="dot"))
fig.add_vline(x=ci_high, line=dict(color="black", dash="dot"))

fig.update_layout(
    title="Jeffreys (arcsine) prior and resulting posterior for a Binomial proportion",
    xaxis_title="p",
    yaxis_title="density",
)
fig.show()

print("Jeffreys posterior mean:", post_jeffreys.mean())
print("Jeffreys 95% credible interval:", (ci_low, ci_high))

Jeffreys posterior mean: 0.16666666666666666
Jeffreys 95% credible interval: (0.04413134197515587, 0.348577710858454)

10.3 Generative modeling (sampling extreme probabilities)

If you sample (p) from an arcsine distribution and then sample data conditional on (p), you generate datasets where the latent probability is often near 0 or 1.

One simple example:

• sample (p \sim \mathrm{Beta}(\tfrac12,\tfrac12))

• then sample a count (K \mid p \sim \mathrm{Binomial}(n,p))

Compared to a uniform prior over (p), this produces more extreme counts (very small or very large (K)).

n = 50
m = 30_000

p_arcsine = stats.arcsine.rvs(size=m, random_state=rng)
p_uniform = rng.random(size=m)

k_arcsine = rng.binomial(n, p_arcsine)
k_uniform = rng.binomial(n, p_uniform)

fig = go.Figure()
fig.add_trace(go.Histogram(x=k_arcsine, histnorm="probability", name="p ~ arcsine", opacity=0.6, nbinsx=n + 1))
fig.add_trace(go.Histogram(x=k_uniform, histnorm="probability", name="p ~ uniform", opacity=0.6, nbinsx=n + 1))
fig.update_layout(
    title=f"Counts from Binomial(n={n}, p) with different priors over p",
    xaxis_title="k (number of successes)",
    yaxis_title="probability",
    barmode="overlay",
)
fig.show()

print("P(k=0)  arcsine vs uniform:", np.mean(k_arcsine == 0), np.mean(k_uniform == 0))
print("P(k=n)  arcsine vs uniform:", np.mean(k_arcsine == n), np.mean(k_uniform == n))

P(k=0)  arcsine vs uniform: 0.08146666666666667 0.0201
P(k=n)  arcsine vs uniform: 0.07906666666666666 0.020866666666666665

11) Pitfalls

• Invalid parameters: always ensure (a<b) (or scale>0).

• Endpoint singularities: the pdf diverges at (a) and (b). This is mathematically fine, but numerically:

  • avoid evaluating the pdf exactly at the endpoints when plotting

  • prefer logpdf for likelihood computations

• Floating-point drift: for the CDF, expressions like ((x-a)/(b-a)) can be slightly below 0 or above 1 due to rounding; clipping prevents sqrt/arcsin from producing NaNs.

• Fitting endpoints: the likelihood can become unbounded as endpoints approach the sample min/max, so naive MLE fitting is ill-posed without additional constraints.

# Numerical illustration: pdf spikes and logpdf is usually safer
x_test = np.array([0.0, 1e-12, 0.5, 1 - 1e-12, 1.0])
print("x:", x_test)
print("pdf:", stats.arcsine.pdf(x_test))
print("logpdf:", stats.arcsine.logpdf(x_test))

x: [0.  0.  0.5 1.  1. ]
pdf: [          inf 318309.886184      0.63662  318313.407023           inf]
logpdf: [      inf 12.670781 -0.451583 12.670792       inf]

12) Summary

• The arcsine distribution on ([0,1]) is (\mathrm{Beta}(\tfrac12,\tfrac12)) with pdf (f(x)=\tfrac{1}{\pi\sqrt{x(1-x)}}).

• It has a U-shaped density: lots of mass near the boundaries.

• Endpoint form on ([a,b]): (f(x\mid a,b)=\tfrac{1}{\pi\sqrt{(x-a)(b-x)}}).

• Mean (=\tfrac{a+b}{2}), variance (=\tfrac{(b-a)^2}{8}); skewness 0, excess kurtosis (-\tfrac{3}{2}).

• Fast exact sampling: (x = a + (b-a)\sin^2(\pi u/2)) with (u\sim\mathrm{Unif}(0,1)).

• In practice, scipy.stats.arcsine gives pdf, cdf, rvs, and fit (with caution near endpoints).